%% Sample Input for Coarse Mesh Diffusion Solver COMEDS2D
% !! This input will correspond to a version which does
%    not use interpolation, generating new responses for
%    each outer iteration !!
% j.roberts 2/23/2010
clear; % clc;
%profile on
name = 'sample_input_v1';

%--------------------------------------------------------------------------
% BASIC PROBLEM INPUT VARIABLES

% number of element types
numtypes = 2;
% number of groups
numg  = 2; 
% number of materials
numm = 2;
% legendre order
order = 3;

%--------------------------------------------------------------------------
%  READ IN DATA -- HERE, JUST DEFINED 
%        D      sigR    nu*sig        scattering
dat = [ 1.500   0.030   0.000        0.000  0.000  
        0.420   0.080   0.135        0.020  0.000
        2.000   0.040   0.000        0.000  0.000
        0.300   0.010   0.000        0.040  0.000 ];    
% absorption cross-section    
abdat = [0.010
         0.080
         0.000
         0.010];
     
%--------------------------------------------------------------------------     
%  PLACE THE ELEMENTS
elements = [ 1  1 2
             1  1 2
             2 2 2];  
% 0.90487828

if (max(max(elements))>numtypes)
    disp('numtypes < largest value of elements')
    return
end

%  Note about numbering: 
%            1 3 5
%            2 4 6 ...etc... bottom to top, since array is sideways

%--------------------------------------------------------------------------
% SET THE BOUNDARY CONDITIONS:  0==vacuum, 1==reflective
%  bc(1) --> left,   bc(2) --> right
%  bc(3) --> bottom, bc(4) --> top
bc(1) = 1;
bc(2) = 1;
bc(3) = 1;
bc(4) = 1;

% modularizing:
%  run geomsetup
%  run solver
%  run output

%--------------------------------------------------------------------------
% SINGLE ELEMENT DISCRETIZATION
% -- this does not necessarily need to be the same for all elements
%    but for now, it's easiest to work with squares
Y=2;
xcm  = [   0     2     4      6  ]*3;
xfm  = [     1      1      1     ]*Y;
ycm  = [   0     2     4      6  ]*3;
yfm  = [     1      1      1     ]*Y;
% compute number in each region
numycm = length(yfm); % number of y regions
numxcm = length(xfm); % number of x regions
mt = zeros(numxcm,numycm,numtypes);

% ELEMENT 1 - material specification
mt(:,:,1)   = [  1 1 1 
                 1 1 1 
                 1 1 1  ]; 
% ELEMENT 2 - material specification
mt(:,:,2)   = [  2 2 2 
                 2 2 2 
                 2 2 2  ]; 
                                      
tol     = [1.e-6 1.e-6]; % tolerance for norms             
maxit   = 15;

swpi = 1;
swnt = 1;
swgm = 0;

input = makeinput( name, numtypes, numg, numm, order, dat, abdat, ...
                   elements, bc, xcm, xfm, ycm, yfm, mt, tol, maxit, ... 
                   swpi, swnt, swgm ) ; 
              
% return  
%-------------------------------------------------------------------------
% RUN SOLVER and OUTPUT EITHER TO FILE OR SCREEN      

tic
[x1,normk1,it,nrm1] = solverPI(input);
fprintf(1,'  final keff = %12.8f\n',x1(end-1));
toc
return
tic
input.maxit=1;
[xseed,tmp,tmp,tmp]=solverPI(input);
[x2,normk2,it2,nrm2] = solverNT(input,xseed);
fprintf(1,'  final keff = %12.8f\n',x2(end-1));
toc

%twoDtwo_sparse_mod
%profile viewer

% possible solver runtime parameters
%     max outers (for iterative)
%     tolerances
%     precondition
%     which solver(s) to use
%

% single element, reflect on 4
%   Y    O     k1                 k2                  %      
%   6    3     1.083028163054769  1.082750633028797   0.025631943081439
%  12    3     1.079839200602866  1.079703240094721   0.012592396048852
% single element, vac on right and top
%   Y    O     k1                 k2                  %      
%   6    3     0.625578053601234  0.625658197825914   -0.012809
%  12    3     0.630559906657183  0.630619579976003   -0.009462
%  two
%  Y    o     k1                 k2                  %     
%  2    3     1.02473269         1.021663636123922  0.300397
%  4    3     1.01584335         1.015112556936627  0.071991
%  8    1     1.01349115         1.013164642941905  0.032226
%twoDtwo_sparse_mod
